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Aims. In this paper, we present a new method to perform numerical simulations of astrophysical MHD flows using the Adaptive Mesh 
f"^ , Refinement framework and Constrained Transport. 

Methods. The algorithm is based on a previous work in which the MUSCL-Hancock scheme was used to evolve the induction equation. In 
this paper, we detail the extension of this scheme to the full MHD equations and discuss its properties. 

Results. Through a series of test problems, we illustrate the performances of this new code using two different MHD Riemann solvers (Lax- 
O ' Friedrich and Roe) and the need of the Adaptive Mesh Refinement capabilities in some cases. Finally, we show its versatility by applying it to 
, two completely different astrophysical situations well studied in the past years: the growth of the magnetorotational instability in the shearing 
CZ2 ■ box and the collapse of magnetized cloud cores. 
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Conclusions. We have implemented a new Godunov scheme to solve the ideal MHD equations in the AMR code RAMSES. We have shown 
that it results in a powerful tool that can be applied to a great variety of astrophysical problems, ranging from galaxies formation in the early 
universe to high resolution studies of molecular cloud collapse in our galaxy. 



' Key words. - MHD - Methods: numerical • 



1. Introduction 

Developing efficient numerical algorithms for the equations of 
magnetohydrodynamics (MHD) is of great astrophysical inter- 
est. Magnetic fields are ubiquitous in a great variety of envi- 
ronments. They are important components of the dynamics in 
such places as the early universe, the interstellar and intergalac- 
tic medium, the environment and interior of stars and the accre- 
tion flow around young stellar objects. 

In the last few decades, finite differences methods have 
been widely used in investigations of a number of astro- 
physical situations in which th e magnetic field is impor - 
tant with su ch codes as ZEUS JStone & Normanl Il992albl) . 
NIRVANA dZiegler & Yorkel Il997h or the Pencil Code 
jBrandenburg & Doblenl2002l) for example. Even though, as 
expected, th e numerical method breaks down in some cir- 
cumstances llFallell2002l) . a considerable amount of progress 
have been made in our understanding of MHD in astrophysics. 
A few attempts have also been made to try to extend the 
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Smoothed Particle Hydrodynamics (SPH) method to M HD 
JPhillips & Monagha jfl985l IPrice & Monaghanll2004albT) . At 
the moment, it is not clear, however, how efficient the resulting 
codes will prove to be in the future. 

In the last few years, several attempts ha ve been ma de to 
try to extend the standard Godunov approach ( lToroll997h . ini- 
tially designed to solve the Euler equations, to MHD. In addi- 
tion to the accurate description of new waves that are peculiar 
to MHD (Alfven waves, the slow and fast modes), one of the 
most dramatic challenge in the development of such schemes 
comes from the solenoidality constraint, which states that the 
divergence of the magnetic field has to vanish everywhere at 
all times. The first algorithms that attempted to solve this prob- 
lem kept the cell centering strategy of the standard Godunov 
approach. T hey used either a "diverge nc e cleaning" step (see 
for example iBrackbill & Barneslll980l or lRvu et alJll998h . or 
various reformulations of the MHD equation s including ad- 
ditional diverge nce-waves (|Powell et alJll999l) or divergence- 
damping terms dDedner et al.ll2002 ) to enforce the solenoidal- 
ity constraint. A novel cell-centered MHD scheme has been re- 
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cently developed bv lCrockett et alJ J2005b that combines most 
of these ideas into one single algorithm. Alternative approach 
used the "staggered" discretisation of the grid commonly 
used in "ZEUS-like" codes along with th e more geometri- 
cal C onstrained Transport (CT) algorithm ( Evans & Hawlev 
1988J). This is for ex ample the case of Balsara & Srjicer 
d 1 999h. iTothl d2000h andlLondrillo & Del ZannaN2000U2004l) . 
iGardiner & StoneH2005al) also explored the possibility of com- 
bining the CT algorithm with the PPM scheme in the new code 
ATHENA. 

Recently, we proposed to extend the well-known MUSCL- 
Hancock algorithm origi nally designed for th e Euler equation 
to the induction equation dTevssier et alJ2006l) . We showed that 
three variants of our scheme have good performances. Two are 
compatible with the Adaptive Mesh Refinement (AMR) algo- 
rithm implemented in RAMSES ( lTevssieri20 Q2). This first part 
was limited to the induction equation, and could only be ap- 
plied to situations where the magnetic field does not affect the 
flow. This is enough, however, to capture the physics of fast dy- 
namos, especially with the help of the AMR. Here we extend 
our approach to the full set of MHD equations and implement 
it in RAMSES. 

The plan of the paper is as follows: in section|3 we present 
the details of the nu merical algorithm. Th e discussion is based 
on our earlier work l lTevssier et al.ll2.Q06h . where the technical 
details of the scheme are presented. In section [3] and [4] we il- 
lustrate the properties of the code on standard ID and 2D test 
problems. In section [5] it is used to study a few 3D flows of 
astrophysical significance: the growth of the magnetorotational 
instability in accretion disks and the collapse of magnetized 
cloud cores. Finally, we summarise the properties of the code 
and highlight future possible developments in section[6] 

2. The numerical method 

2.1. Equations and notations 

The equations we seek to solve are the usual MHD equations. 
When written in conservative form, they read: 



^+V-(pv) = 0, 

ot 

^L+V-{pvv-BB) + VP lot = 0, 
Ot 

^+V- [(E + P tot )v-B(B-v)] = 0, 
ot 

dB 

— + V-(vB - Bv) = . 

ot 



(1) 

(2) 
(3) 
(4) 



Here, p is the fluid density, v its velocity and B is the magnetic 
field. P to t stands for the total pressure, the sum of the thermal 
pressure P and the magnetic pressure: 



P,n, =P + 



BB 

~2~ ' 



(5) 



and E is the total energy of the fluid 



E = e 



vv BB 



where e denotes the internal energy. Unless otherwise stated, 
we will assume throughout this paper that the equation of state 
is that of a perfect gas, in which case P = (y — l)e. 

As discussed in the introduction, this set of equations has 
to be completed by the solenoidal constraint, to be satisfied at 
all times: 

VB = (7) 

As in iTevssier et al. (2006), we will use t hroughout this pa- 
per the CT scheme fevans & Hawlevl ll988) to ensure that this 
condition is fulfilled to machine-roundoff precision. It simply 
consists in writing the induction equation in integral form: 



E-dl. 



(8) 



where E is the electric field defined by the relation E = vxB. 
While all the hydrodynamic variables (density, velocities, total 
energy) are located at cell centers, this approach requires the 
magnetic field components to lie on the cell faces. The grid 
structure that results is described in the following section. 

2.2. The staggered mesh 

In the followings, we describe our scheme using 3 dimensional 
coordinates x, y and z. The physical variables are discretized 
on a standard 3D Cartesian grid. The center of each cell is lo- 
cated at the position (xj,yj,Zk)- In a given cell, faces normal to 
the x-direction have coordinates x = jc, ± i/2 and cover a surface 
element defined by yj-i/i < y < X/+1/2 and Zk-i/2 < z < Zk+i/2- 
The coordinates of the other faces, normal to the y and z direc- 
tion, can be similarly defined. 

As for the Euler equations, the hydrodynamical variables 
(density, momentum, energy) are volume-averaged over a cell 
and the discretized values are defined at the cell center. For 
example: 



1 



Pi,j,k 



AxAyAz 



rxt+i/2 ryj+i/2 rzk+m 

p(x',y',z')dx'dy'dz' (9) 

JXi-l/2 Jyj-1/2 ^Zk-i/2 



Because of the staggered mesh representation, magnetic fields 
components are surface-averaged over the cell face to give: 

B x ,i-u2J,k=j-j- B x {x i ^ l2 ,y',z')dy'dz' (10) 

Jyj-l/2 •JZt-1/2 

Here Ax, Ay and Az stand for the Cartesian mesh size in each 
direction. 



2.3. The Euler system 

As outlined in lLondrillo & Del Zannal l l2000l) . the system of 
MHD equations written in section 12.11 can be broken in two 
sub-systems. The first involved the time evolution of the cell- 
centered, volume-averaged variables and is reminiscent of the 
standard Euler equations, which includes mass, momentum 
and energy conservation. This set of equations, quite naturally 
called the "Euler system", can be written in vectorial form 



dU dF dG dH „ 



(11) 
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where 

U = (p,pv x ,pv y ,pv-,E) T 

and the flux function F is given by 



(12) 



variables, we need to define volume-averaged magnetic field 
components. We use for that purpose the average of their cor- 
responding face-centered components 



F = 



PV X 

• 2 - D - Bl 



B ",i,j,k ~ ? { B ",i-l/2,j,k + B x,i+l/2J,k) ' 



(16) 



pvj + P tc 

pV X Vy 

pv x v z 



B x B y 
BJB, 



(E + P tot )v x - B x (B-v) 



(13) 



The expression for G for H are completely symmetric. 
Integrating in space over a cell and in time between f and f 1 , 
equation (II It writes: 



U" + l -I/"., 

i.J,k i,j,k 
At 



-iK+1/2 



7 n+l/2 



i+l/2,7,* x i-l/2,j,k 



Ax 



and likewise for B"..^ and B". We finally augment the vec- 
tor W" with these 3 new components. We use in our cur- 
rent implementation two standard slope limiters (used in many 
fluid dynamics codes), the MinMod slop e and the MonCen 
(Monotonized Central) slope JToroll 19971) . The MinMod lim- 
iter is more diffusive than the MonCen limiter, so we use the 
latter for most applications. On the other hand, the MinMod 
limiter is known to ensure the positivity of the solution in mul- 
tiple space dimensions. In difficult cases, we therefore switch 
to the MinMod slope limiter. 



G 



n+l/2 



Qn+l/2 



ij+1/2,* "/j- 1/2,* 



Ay 

jjii+1/2 _ zjn+l/2 
+ "i,/,*+l/2 "/,/,*- 1/2 _ 

Az 



(14) 2.3.2. The Euler predictor step 



where superscripts n and n+ 1 refer respectively to time coor- 
dinates f and f" +1 . U" and t/" +1 are the volume-averaged vari- 
ables at time t" and f" +1 . The time- and surface-averaged fluxes 
are defined by 



The standard MUSCL methodology can be applied to the Euler 
sub-system, using the previously defined cell-centered vector 
W. The solution is advanced in time up to t n+l11 using a Taylor 
expansion of the Euler system in non-conservative form based 
on the previously computed TVD slopes 



pn+l/2 
1+1/2,./,* 

1 

AtAyAz 



W «+W _ W n _ A n 

n i,j,k ~ n i,j,k A *i,jjc\ Q x 



At/2 



,.j.k 



ldW\" (dW\" 



Similar time and surface-average quanti- 

«+l/2 , TT n+ 1/2 

i,; + i/2,* and H i,j,k-H2 appearing 



where At = f" +1 - f 
ties are written for the fluxes G 
in equation Jl 5I >. 

In this paper, w e intend to extend the well-kn own MUSCL- 
Hancock scheme Jvan Lee J Il977t iTorol Il997l) to the equa- 
tions of ideal MHD. When applied to the Euler equations, 
this method performs the conservative update of the volume- 
average variables U in two steps: a predictor step and a correc- 
tor step. In the former, the vector U is computed at the half time 



step f 



H+l/2 



= t" 



At 12 using a Taylor expansion of the under- 



lying hyperbolic system. It is also spatially reconstructed from 
the cell center to the cell faces using a piecewise linear recon- 
struction based on TVD slope limiters. In the predictor step, 
the fluxes appearing in equation dl5> are evaluated by solving 
a ID Riemann problem between the two (left and right) recon- 
structed states at each cell interface. 

2.3.1 . Cell-centered TVD slopes 

The first step in the MUSCL approach is the computation of 
finite-difference approximation of the spatial derivatives of all 
cell-centered quantities. As usually done in higher order finite 
volume schemes, spatial derivatives are approximated using 
slope limiters, in order to obtain positivity preserving, non os- 
cillatory solutions. Except for the final conservative update, we 
always use the primitive variables W — (p, v x , v y , v z , P) T in all 
intermediate calculations. In addition to these 5 cell-centered 



where the matrix A x (resp.A^, and A z ) is the Jacobian matrix 
of the flux function in the x (resp. y and z) direction, evalu- 
ated using the cell-averaged state W"j k . At this stage, we have 
cell-centered predicted states at time f +x l 2 for the 5 Euler vari- 
ables W"^ 2 = (p,v x ,Vy,v z ,P) T . Face-centered predicted val- 
ues for the magnetic field components are also computed using 
the method described in details in section l2~4l We compute the 
predicted cell-centered components of the magnetic field using 
the average of their corresponding face-centered values. We fi- 
nally augment the vector W"^ l J 2 with these 3 new cell-centered 
predicted variables. 

2.3.3. The Euler corrector step with 1 D Riemann 
solvers 

Using the TVD slopes computed at time f, we reconstruct the 
primitive variables at each cell-interface, except for the longi- 
tudinal magnetic field component, since its predicted value has 
been already computed at the correct location (see Sect. I2.4t . 
For example, at the two interfaces perpendicular to the x-axis, 
we obtain the two following reconstructed states 



W 



«+l/2,L 
i+l/2J,* 



dx )i,j,k 



w 



n+l/2,« 
i-U2,j,k 



w 



n+l/2 

hJJt 



9W\" . . 
IT Ax/2 

8X li.Lk 



(18) 



(19) 
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These states will be used as input states for ID Riemann prob- 
lems perpendicular to each interfaces. Note that for these ID 
MHD Riemann problems, left and right states are defined us- 
ing only 7 variables, namely the 5 Euler variables and the 2 
magnetic field transverse components, thus the name "seven 
waves Riemann solvers". As far as the Riemann solver is con- 
cerned, the longitudinal component of the magnetic field is as- 
sumed to be constant in time and space, in order to enforce the 
solenoidality constraint in one space dimension. This constant 
value is taken equal to the predicted value at time f" +1 / 2 , namely 



B" +1 J. 2 n .. (resp. Z?" + V z . „ , and B"TV* , ,,) for the interface per- 

x,/+l/2j,A- v " y,i,j+l/2,k z,i,J,k+i/2' " 

pendicular to the x (resp. y and z) axis. The output of these ID 
Riemann solvers are the time and surface averaged fluxes at the 
same interface F^ jJc (resp. G${^ and fl££f 1/a ). In our 
current implementation, we use two different Riemann solvers, 
namely a simple, local Lax Friedrich (LLF) solver, for which 
the flux is given by 

Fllf(W l , W r ) = 1 -{F L + F R ) - i max |^| (U R - U L ) (20) 

2 2 a=l,7 



yi+1/2 



,n+l/2 



and the MHD Roe solver described in ICargo&Gallicel(ll997l) 
and developped by iGardiner & Stond 112005 ah . for which the 
flux can be written as 



F Roe (W L ,W R ) = l -(F L + F R ) 



- 5 2 Ra\A a \L a ■ (U R - U L ) 



(21) 



<r=l,7 



where Ui and Ug are the conservative state on each sides of 
the interface, Fi and Fg the associated fluxes, R a and L a are 
respectively the right and left eigenmatrices of the Roe matrix 
and A a its eigenvalues (wave speeds). 

2.4. The induction system 

To form the full set of MHD equations, equation il Q has to be 
completed by the induction equation, called here the induction 
sub-system. As for the Euler system, the induction equation 
can be written in conservative form by a straightforward inte- 
gration in space-time 



«+ 1 



B" + 



- B" . 



7 n+l/2 



-E 



n+1/2 



x,i-l/2J,k x,!-l/2J,jfc e,/-l/2,/+l/2,A' ^z,i-l/2J-\/2,k 



At 



Ay 



(22) 



E-n+l/2 _ pn 

y,i-l/2,j,k+V2 £ 'y./-l/2.>,A-l/2 _ 

Az 



7 n+l/2 



with similar expressions for B" +l and B" +l . Here conventions 
are similar to the ones used in section l2~3l above. except that one 
defines now a time- and edge-averaged electromotive force 
(EMF) as 



^i+l/2 



~z,;-i/2,./~i/2,* 

1 /"*'''+! rZk+l/2 

j^ z j J E z (xi- m ,yj-ii2,z',t')dz'df , 



(23) 



and similar expressions can be derived for E" +l ^ 2 and E" +1 ^ 2 . 
As for the Euler system, the numerical evaluation or the EMF 



proceeds in two steps: a predictor step followed by a correc- 
tor step. The MUSCL methodology can be extended to the in- 
duction system and t his extension was extensively discussed in 
iTevssier et alJ l EoOrjl) . We recall here only the basic ingredients. 

2.4.1 . Face-centered TVD slopes 

In order to obtain a second-order accurate and non-oscillatory 
solution, we need to use spatial reconstruction of the magnetic 
field components based on TVD slope limiters. The main dif- 
ference arises because of the finite-surface representation of 
the magnetic field. Indeed, we need a piecewise linear repre- 
sentation of B A within the y-z plane. For that purpose, we use 
the same TVD slopes (MinMod and MonCen) as above, using 
the face-averaged value of the 3 magnetic field components at 
time f" . For B x , we need to compute only the 2 transverse slopes 
dB Jdy and dB x jdz. A similar property holds for B v and B-. 

2.4.2. The induction predictor step 

Various methods to perform the predict or step for the induc- 
tion system were recently explored by ITevssier et alJ 
for the kinematic case. These methods were referred to as 
Runge-Kutta, U-MUSCL and C-MUSCL. The extension of 
the first two to the full set of MHD equations, while possi- 
ble, is computationally expensive because they require to solve 
one (U-MUSCL) or two (Runge-Kutta) Riemann problems in 
the predictive step. Moreover, the large stencil of the Runge- 
Kutta scheme is not compatible with the compact stencil re- 
quired by our tree-based AMR implementation. For these var- 
ious reasons, we decide to use only the C-MUSCL scheme 
in our current MHD application. It combines the nice prop- 
erties of being computationnally efficient and compatible with 
the AMR requirements. The price to pay is a reduced stabil- 
ity range for the time step, since the Courant factor has to be 
less than 2/( a/2 + 1 ), instead of 1 for the other schemes (see 
ITevssier et alJ200rl for details) . 

The purpose of the predictive step is to advance the solution 
between t" and f +x l 2 using the CT algorithm. For that, EMF 
need to be spatially interpolated on cell edges at time t". The 
idea of C-MUSCL is to do it by simple arithmetic averages of 
the magnetic field and velocity components. For example, the 
EMF E^ f _y2 j-1/2 k ^ s calculated by: 



E l>-i/2j-i/2,k ~ v* B y ~ Vy B *> 
with 

V ~y = 5<< ; . ; i-'v, I.U + <,; U + l v, I./ 

^ x ~ x,i-\l2,j,k + B x,i-l/2,j-lj) ' 
1 



By - y( B y,iJ-\/2,k + K,i-\J-\/2,k) ■ 



(24) 

(25) 
(26) 
(27) 
(28) 



This spatial reconstruction is second order in space, altough 
TVD slopes have not been used at that time. The EMF 

£ ",i,/-l/2,*-l/2> E 'y\i-l/2J,k-l/2 and E li-V2,j-l/2,k are then USed t0 
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update the solution between t„ and f n+ i/2 using the CT algo- 
rithm (see Eq. I22> . Because only one EMF is calculated per 
cell edge, the predicted face-centered magnetic field (B" +l ^ 2 , 
B" y +l ^ 2 and B" +x ^ 2 ) satisfies the solenoidality constraint ex- 
actly. The properties of C-MUSCL and their comparison with 
the U-MUSCL and Runge-Kutta schemes are described in 
iTevssier et all d2006l) . It was found that C-MUSCL behaves 
essentially similarly to these two other schemes, with a lower 
computational cost and a slightly stronger Courant condition. 



2.4.3. The induction corrector step with 2D Riemann 
solvers 



As described above, after the predictor step, we have obtained 
the 5 cell-centered Euler variables W"^ 2 = (p, v x , v y , v,, P) T 

and the 3 face-centered magnetic field components B" +l ^ 2 , 
fi" +1 ^ 2 and B" +l ^ 2 . Using equation at time f" +1 ^ 2 , we have 
also obtained the 3 predicted cell-centered components of the 
magnetic field. We finally augment the vector W"^ 2 with these 
3 new cell-centered predicted variables. 

Following again the MUSCL methodology, we now need 
to reconstruct complete MHD states at each cell-edge, in order 
to compute the EMF for the final conservative update of the 
magnetic field components. This reconstruction will produce 4 
different states in the 4 cells adjacent to the edge. For obvious 
reasons, these 4 states, separated by 4 boundaries (labelled N 
for North, S for South, W for West and E for East) will be 
labelled in clockwise order by NE, SE, SW and NW. 

We first reconstruct the cell-centered state to the cell edges 
using 



■n+l/2,NE 



w 



n+1/2 
ij,k 



3WY Ax_(dW\' Ay (29) 



8y 



and similar relations that defines W"^j 2 'j_ E /2k , W"^j 2 'j^ /2k 
and W^lip'j-ipk' S mce tne 2 longitudinal magnetic field com- 
ponents (B" +1/2 and B" +1/2 ) are already defined at the 4 adjacent 
interfaces, we only need the cell-centered transverse compo- 
nent B" +1/2 in the above reconstruction. B" +l/2 and B" +1/2 are 
reconstructed using face-centered TVD slopes as 



B 



n+l/2,S 
x,i-l/2,j-\/2,k 



= B 



,n+l/2 

x,i-\/2,j-\,k ' 



B 



,n+l/2,W 
y,i-l/2,j-U2,k 



B' 



,«+l/2 

y,;-l,./~l/2,* ' 



dBX Ay 

°y h-i/2,j-i,k 2 
8B y \ n Ax 



d x 'i-l,j-l/2,k 2 



(30) 



(31) 



and similar relations defining B" + ,,\ . , ,~ , and B" + 1.1 . , ,, , . 
The four corner states, with 6 variables each, and the 4 longi- 
tudinal magnetic field components entirely define a 2D MHD 
Riemann prob lem, which satisfies the soleno idality constraint 
in a 2D sense. lLondrillo & Del Zannal feOOOl) have shown that 
the EMF entering in the final Contrained Transport update 
should be obtained as the solution of this 2D Riemann problem, 
in order to obtain a stable numerical solution, with a proper up- 
winding of all MHD waves. While a very simple exact solu- 
tion exists in the kinetic case ( ITevssier et alJl200rj|) . designing 
2D Riemann solvers for the full set of MHD equations is an 



ambitious task that is beyond th e scope of this paper. An ap- 
proximate sol ution, proposed bv lBalsara & Spice J ll 19991) and 
IZiegleiU2004h . is based on averaging the flux given by the four 
adjacent ID Riemann problems. The solution E 2D of the 2D 
Riemann problem writes in that case: 



E 2D (W NE , W SE , W sw , W NW ) 



-(E W (W NW ,W NE ) + E 1D (W SW ,W SE ) 



+El D (W sw ,W NW ) + E iU (W^,W"*)) 



(32) 



where the quantities El D stands for the solution of the ID 
Riemann problems defined by the four states. This solution re- 
lies on 4 Riemann solvers per cell ed ges and turns out to be 
quite expensi ve. Following the ideas of lLondrillo & Del Zannal 
(200(1 12004) . we exploit the fact that in our current imple- 
mentation, we use only linear Riemann solvers (namely Lax- 
Friedrich and Roe). In this case, the flux can be written as in 
equation!^ If we now use 2 Roe matrices, one for each direc- 
tion, instead of 4, the EMF function can be written as 

E 2D (W NE , W SE , W SW ,W NW ) = 

-(E Z (W NE ) + E Z (W SE ) + E Z (W SW ) + E-(W NW )) 



-\YjKK\K ■ (Ue-u w ) 

<i=1,7 

+\ Yj ■ (Un-Us) 



(33) 



a=lJ 



where Ue, Uw, Un an d Us are averaged conservative vari- 
ables defined at the interfaces of the four ID Riemann prob- 
lems. One intersesting property in the above expression is the 
explicit contribution of the 2 diffusive terms coming from the 
2 Roe matrices. For non-linear Riemann solvers, such as HLL, 
it is preferable to use the 4 Riemann solv ers approximation 
(lLondrillo & Del Zannal2004HZiegleJ2004 . 

2.5. The AMR scheme 

T he AMR algorith m used in RAMSES is descr ibed in details 
in Tevssie rl ll2002l) and its extension to MHD in ITevssier et alJ 
yOOfJ). We briefly recall the main features here. It is a tree- 
based AMR code originally designed for cosmologic al applica- 
tions. The data structure is a "Fully Threaded Tree' ' JKhokhlovl 
119981) . The grid is divided into groups of 8 cells, called "octs", 
that share the same parent cell. Each oct has access to its par- 
ent cell address in memory, but also to neighboring parent cells. 
When a cell is refined, it is called a "split" cell, while in the op- 
posite case, it is called a "leaf" cell. The computational domain 
is always defined as the unit cube, which corresponds in our ter- 
minology to the first level of refinement in the hierarchy C — 1 . 
The grid is then recursively refined up to the minimum level of 
refinement {,„,„, in order to build the coarse grid. This coarse 
grid is the base Cartesian grid, covering the whole computa- 
tional domain, from which adaptive refinement can proceed. 
This base grid is eventually refined further up to some maxi- 
mum level of refinement { max , according to some user defined 
refinement criterion. When ( max = £ m j n , the computational grid 
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is a traditional Cartesian grid, for which the previous scheme 
apply without any modification. When refined cells are created, 
however, some issues specific to AMR must be addressed. 

2.5.1. Divergence-free Prolongation Operator 

When a cell is refined, eight new cells (i.e. a new "oct") are 
created for which new cell-centered variables and new mag- 
netic field components are needed. This operation is usually 
referred to as the "prolongation operator". The traditional ap- 
proach relies on a conservative interpolation of the 5 cell- 
centered conservative variables U = (p,pv x ,pv v ,pv z , E) T . For 
the face-centered variables, each of the six faces of the parent 
cell are split into 4 new fine faces. Three new faces, at the cen- 
ter of the parent cell, are also split into four new children faces. 
The resulting magnetic field components, fine or coarse, need 
to satisfy the divergence-free constraint in integral form. 

This critical st ep has been solved by iBalsaral J200ll) and 
lT6fh&RoeH2002l) in the CT framework. We recommend both 
of these articles for a detailed description of the method. The 
idea is to used slope limiters to interpolate the magnetic field 
component inside each parent face, in a flux-conserving way, 
and then to use a 3D reconstruction, which is divergence-free in 
a local sense inside the whole cell volume, in order to compute 
the new magnetic field components for each central children 
faces. In our case, the same slope limiters as in the Godunov 
scheme (MinMod or MonCen) can be used. 

This prolongation operator is used to estimate the mag- 
netic field in newly refined cells, but also to define a tempo- 
rary "buffer zone", two "ghost cells" wide, that set the proper 
boundary for fine cells at a coarse-fine level boundary. This is 
the main reason why compact stencils are needed for the un- 
derlying Godunov scheme. 

2.5.2. Magnetic Flux Corrections 

The other important step is to define the reverse operation, 
when a split cell is de-refined, and becomes a leaf cell again. 
This operation is usually called the Restriction Operator in 
the multigrid terminology. The solenoidality constraint needs 
again to be satisfied, which translates into conserving the mag- 
netic flux. The magnetic field component in the coarse face 
is just the arithmetic average of the 4 fine face values. This 
is re miniscent of the "flux correction step" for the Euler sys- 
tem (IBerger & Qligerlll984l IBerger & Colellalll989t iTevssierl 

2.5.3. EMF Corrections 

The "EMF correction step" is more specific to the induction 
equation. For a coarse face which is adjacent, in any direction, 
to a refined face, the coarse EMF in the conservative update of 
the solution needs to be replaced by the arithmetic average of 
the two fine EMF vectors. This guarantees that the magnetic 
field remains divergence-free, even at coarse-fine boundaries. 
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Fig. 1. Amplitude (upper panel) and phase velocity (lower 
panel) of the circularly polarized Alfven wave as a function of 
time. The full lines display result using a Roe solver whereas 
the dotted lines show results obtained with a Lax-Friedrich 
solver. The resolutions are from top to bottom, 100, 30 and 10 
grid points per wavelength 



3. Numerical tests in 1D 

3.1. Non-linear Alfven wave test 

The first test we present is the propagation of non-linear cir- 
cularly polarised Alfven waves. Such waves, which are exact 
solution of the MHD equations, propagate in a gas of uniform 
density, po and along a uniform magnetic field, Bq z . They are 
given by: B x - B ± cos(u>t - kz), B y = B ± sin((L>t - kz), V x = 
V ± cos(a>t - kz), V y = V±sin(o>t - kz) where o)/k = Bq z / j4-npo 
andB ± /V ± = y/4np . 

We have simulated the propagation of these waves on a uni- 
form grid, for Bq z = V47T, B ± = an d k = 0.1. This leads 
to a wave period equal to 0. 1 . The agreement between the an- 
alytical and numerical solutions depends on the numerical res- 
olution. Fig.[2displays the wave amplitude (upper panel) and 
phase velocity (lower panel) as a function of time for the Roe 
and Lax-Friedrich solvers and different resolutions, namely 10, 
30 and 100 points per wavelength. With 10 grid points per 
wavelength, the amplitude quickly decays because of numer- 
ical dissipation and in about 5 wave periods, the amplitude of 
the waves is only 40% of its initial value. With 30 and 100 grid 
points per wavelength the agreement is much better and almost 
no decay has occurred even after 50 wave periods in the lat- 
ter case. In the lower panel, the wave velocity is also seen to 
agree better and better with its theoretical value, v„. = 1, as the 
resolution is increased (note that for 100 grid points, the wave 
velocity obtained with the Roe solver is not represented as it is 
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Fig. 2. Solution to the MHD shock tube showing the density as 
a function of x, obtained with RAMSES when a — n. The code 
uses 400 grid points in that case and the Roe Riemann solver. 



indistinguishable from the Lax-Friedrich results). As expected 
the Roe solver leads to slightly better results in both cases than 
the Lax-Friedrich solver. 



3.2. MHD shock tube test 

An interesting application of the AMR scheme is the study 
of the development of compound waves in shock tube calcu- 
lations. It has b een analysed with finite volume schemes by 
lTorrilhonl ( E004l) . through the analysis of the MHD shock tube 
whose initial state is: 



W L = (1,0,0,0, 1,1,1, of 

W R = (0.2,0,0,0,0.2, 1, cos a, sin af 



(34) 
(35) 



where W = (p, v x , v y , v z , P, B x , B y , B Z ) T . When a — n, there 
are two solutions to the Riemann problem: the first is regu- 
lar, which means that it contains only shocks and contact dis- 
continuities. The second, however, features a comp ound wave 
which is a composition of an Alfven and a slow wave. lTorrilhonl 
J2004 showed that finite volume schemes converge toward 
the second. When a is different from n, the solution is regu- 
lar and sh ould only contain shocks and contact discontinuities. 
However. ITorrilhor] J2004I) found that finite volume codes still 
tend to exhibit the compound wave for low and moderate reso- 
lutions. The solution converges toward the regular solution only 
when very large resolutions are used. 

Here, we use RAMSES to illustrate how the AMR scheme 
can help to solve this problem. Figure |2] shows the density 
vs. position when a — n at time t = 0.4. The grid is com- 
posed of 400 cells evenly distributed between -1 and 1.5. 
The AMR is switched off in this first run. The compound 
wave is clearly visible at i ^ -0.25. The whole solution 
looks identical to previous results published in the literature 
with similar codes jRvu & Jonesll995UCargo & Gallicel 19971 
lLondrillo & Del Zannal2000l) . 

Taking now a = 3 and computing the solution of the MHD 
shock tube on a uniform grid, we also found that the compound 
wave remains for low resolution as described by ITorrilhor] 



Fig. 3. Zoom on the region in which the compound wave devel- 
ops when a = 3. Black curves are obtained on a uniform grid. 
From top to bottom (at x ~ -0.24), they correspond to 800, 
1200, 1600, 2000, 3000, 5000, 10000 and 20000 grid points re- 
spectively. The red curve was calculated after switching on the 
AMR scheme and using a similar CPU time as for the 20000 
grid points curve. Its maximum resolution is equivalent to us- 
ing about 10 6 cells on a uniform grid and show a dramatic im- 
provement toward the regular solution. 
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Fig. 4. Complete solution of the MHD shock tube when a = 
3, with the AMR scheme turned on. The solid line is a plot 
of the density as a function of position while the dashed line 
(whose scale is given on the right axis) illustrates the level of 
refinement the code uses for each cell. 



d2004 . This is illustrated in figure which is a zoom on the 
structure of the solution in the neighbourhood of the compound 
wave. The black lines are computed on a uniform grid and cor- 
respond to increasing resolution. Namely, from top to bottom 
(at x ~ -0.24), the number of cells are 800, 1200, 1600, 2000, 
3000, 5000, 10000 and 20000. The red line shows the result 
of the same model computed using the AMR scheme with a 
refinement strategy based on the magnitude of the gradient of 
all 7 flow variables. The finest resolution in this run is equiva- 
lent to having 10 6 cells on a uniform grid. We found the result 
of this model to be almost indistinguishable from the regular 
solution (remember that this is the ONLY physical solution to 
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Fig. 5. Time history of the magnetic energy in the magnetic 
loop advection test. The different curves are obtained using the 
Roe Riemann solver (solid line) and the Lax-Friedrich solver 
(dashed line). 



Fig. 6. Snapshot of the density at time t — 0.5 resulting from the 
Orszag-Tang test. The grid is uniform and composed of 512 2 
cells. The result is computed using the Roe Riemann solver 
and the C-MUSCL predictive step. 



this problem). Interestingly, this is not the case for the uniform 
runs. Even though the compound wave is seen to gradually dis- 
appear as the resolution is increased, features departing from 
the correct solution are still observed even when 20000 grid 
zones are used. This illustrates the extremely large resolution 
needed to accurately calculate the solution of this problem and 
shows the interest of using AMR. Indeed, the AMR run used 
only 10000 cells for the same equivalent resolution as 10 6 grid 
cells, which corresponds to a gain of about 2500 in CPU time. 

The complete solution of the AMR model is shown on fig- 
ure @] As is figure |2 the solid line represents the density as a 
function of position. The dashed line shows the corresponding 
refinement level. It scale is indicated on the right axis. As ex- 
pected, the grid is highly refined at the location of the shocks 
and contact discontinuities. It is important to understand that 
the AMR is not just a fancy tool for this test, but is actually es- 
sential to solve it properly. One might indeed think that increas- 
ing the order of the numerical scheme wo uld help to converge 
to the regular solution at lower re solution . iTorrilhon & Balsaral 
actually showed that the improvement when using third 
or fourth order WENO schemes is small. This is because the ac- 
curacy of any such schemes breaks down to first order close to 
discontinuities, which is precisely where the compound wave 
lies. 

4. Numerical tests in 2D 

4.1. Advection of a magnetic loop 

As a first 2D test, we now consider the simple advec- 
tion of a magneti c loop t hat has recently been proposed by 
iGardiner & Stonel ( l2005ah . It simply consists in the evolu- 
tion of a weak magnetic field in an initially uniform veloc- 
ity field. Since the thermal pressure is much larger than the 
magnetic pressure, the magnetic field can be considered as a 
passive tracer advected in a ti me independent flow. The in itial 
setup is exactly the same as in IGardiner & Stonel J2005al) and 
iTevssier et alJ ( l20o3) . The velocities are set up to v x = 2, v y = 1 



and v- = 0. The initial magnetic field is such that B z = 0, while 
the components B x and B y are defined using the z-component 
of the potential vector A (B = VxA): 



A. = 



A (R-r) forr<R, 







otherwise , 



(36) 



with Aq = 10~ 3 , R = 0.3 and r = ^x 2 + y 2 . The computational 
domain is defined as -1 < x < 1 and -0.5 < y < 0.5. There 
are 128 cells in the first direction and 64 in the second. The 
solution is evolved between t — and t — 2 and we analyzed 
the results obtained by our scheme using the MonCen slope 
limiter, comparing explicitly the Roe solver to the local Lax- 
Friedrich solver. 

A simple way to evaluate the efficiency of the scheme 
is to compare the time history of the magnetic energy E m . 
This is done in figure [5] where E m is represented as a func- 
tion of time for the Roe solver (solid line) and for the Lax- 
Friedrich solver (dashed line). We first note that the results ar e 
very similar to those published bv IGardiner & 810*1161 d2005al) . 
This demonstrates that, using only TVD linear reconstruction, 
our scheme provide s comparable accuracy to the piecewise 
parabolic scheme of IGardiner & Stone] J2005al) . As expected, 
the Lax-Friedrich solver is more dissipative, while the results 
obtained using the Roe solver exact ly reproduce the res ults 
obtained in the pure advection case JTevssier et al.ll2006l) . to 
which we refer the reader for more details, especially regard- 
ing the AMR scheme. 

4.2. The Orszag-Tang vortex 

One of the most well-known 2D MHD test is the Orszag-Tang 
test. The initial condition of the flow create a vortex that is un- 
stable and quickly breaks down into turbulence. Although no 
analytical solution is known for this test, it has been so widely 
documented in the literature that it is now very useful as a first 
2D benchmark for a code. 



Fromang et al.: A High Order Godunov Scheme for Astrophysical MHD 



9 



The initial state is denned as 
P = yPo, 

V = (- sin 2ny, sin 2nx) , 
B = (-Bo sin27rx, Basin Airy) . 
The different parameters are defined by 



(37) 
(38) 
(39) 



grid. We use a refinement strategy based on the gradient of the 
thermal pressure: 



Yin 



and B, 



1 



o 



V47T 



The grid extends from to 1 in both directions and we use 
512 2 cells uniformly distributed over the computational do- 
main. The solution is evolved between t — and t = 0.5 us- 
ing the Roe solver and the MonCen slope limiter. The den- 
sity distribution in the {x,y) plane at the end of the simu- 
lation is shown on figure [6] The agreement between our re- 
sult and previous published w ork is excellent llRvu et all 1 998: 
ILondrillo & Del Zannal2000l) . The complex pattern of interact- 
ing waves is perfectly recovered. 

4.3. The current sheet 

iGardiner & Stone) d2005al) recently described a 2D test that fol- 
lows the time evolution of a current sheet created by a magnetic 
field discontinuity. Reconnection occurs at the location of the 
discontinuity. Because no explicit dissipation is included in the 
code, the entire evolution is driven by the numerical resisitivity 
of the scheme, and, as such, is sensitive to every details of the 
algorithm. The initial setup is described in the followings. 

The computational domain lies in the domain < x < 2 
and < y < 2 and is divided in 256 uniform cells in each 
directions. At time t — 0, density and pressure are uniform: 
p — 1 and P = 0.1. The magnetic field components vanish in 
the x and z direction and B y is defined by 



By = 



-1 if \x- 1| < 0.5, 
+ 1 otherwise . 



(40) 



Similarly, \\, — v- = and v x - \'o sin (ny). 

We first present our numerical results using the Roe 
solver with the MonCen slope limiter. They are represented 
on figure Q The magnetic field lines are plotted at times 
t - 0, 0.5, 1, 15, 2, 2 .5, 3, 3.5 and 4. As reported by 
IGardiner & Stonel J2005al) . magnetic islands form, grow and 
eventually merge with each other. At the end of the simulation, 
four islands are clearly visible at the location of each discon- 
ti nuity. A direct comparison between our results and figure 12 
of IGardiner & Stonel J2005al) shows that both codes agree al- 
most perfectly up to time t = 2.5. On the other hand, at later 
time, no strong evolution is observed in our case, while for the 
ATHENA code, the flow symmetry is broken and the two is - 
lands merge into a single large one (IGardiner & Stonel2005al) . 
As discussed above, this difference is an indication that both 
codes, although very similar, have different dissipative proper- 
ties. 

The next step is to test our AMR scheme in 2D: we perform 
the same exact simulation, except that now we use a base grid 
of n x = n y = 32, which corresponds to l min = 5, with 3 addi- 
tional levels of refinement, so that C max = 8. The formal resolu- 
tion is thus equivalent to the first test with a regular Cartesian 



max (\A X P\ , |Aj,P|) 



> 0.05 , 



(41) 



associated with a similar criterion based on B y . In Figure [8] 
we show the magnetic field lines at time t = 2 obtained with 
the Roe solver (upper right plot): it is indistinguishable from 
the previous result obtained with a Cartesian grid. For sake of 
comparison, we have also compared the result obtained at t — 2 
with the Lax-Friedrich solver (upper left plot): it is now com- 
pletely different. Due to the increased numerical diffusivity of 
the Lax-Friedrich solver, the tiny magnetic islands, that move 
up and down from the center of the image, have not yet merged 
to their final position. Moreover, with the Roe solver, we also 
obtain a static magnetic island at the position of the flow stag- 
nation point. This static island is absent from the Lax-Friedrich 
solution. As anticipated, this test is highly discriminant of the 
dissipative properties of numerical schemes. 

In order to study the convergence of each solution, we take 
advantage of the speed-up provided by the AMR grid to per- 
form additional high-resolution simulation. Using the same re- 
finement strategy, we now set { max = 10, so that the formal 
resolution is now n x = n y = 1024. We present in Figure [9] 
the resulting AMR grid, together with a grey-scale image of 
the thermal pressure. We see that AMR cells are optimally dis- 
tributed in order to properly sample the current sheet, as well 
as sharp MHD waves propagating perpendicular to the current 
sheet. The results obtained at t — 2 using both Riemann solvers 
are shown at the bottom plots of Figure [8] The Roe solution 
has not changed when compared to its low-resolution coun- 
terpart. This demonstrates that the lower numerical dissipation 
of the Roe solver allows a faster convergence of the numerical 
solution. Interestingly, the Lax-Friedrich high-resolution solu- 
tion has also converged toward the same solution, except for the 
static magnetic island at the flow stagnation point. This demon- 
strates that using AMR can provide fast convergence towards 
the true solution, even with a rather dissipative scheme. On the 
other hand, the static island in the center of the flow seems to 
be highly sensitive to the details of the scheme. As opposed 
to Lax-Friedrich, the Roe solver has the interesting property 
that for a static velocity field, the numerical dissipation van- 
ishes exactly. The difference between the 2 solvers is there- 
fore maximized at the stagnation point, where both schemes 
are converging towards 2 different solutions, even at our high- 
est resolution. This peculiar behavior is due to the fact that this 
reconnection problem is performed without any physical resis- 
tivity. It should be therefore considered only as an interesting 
numerical test, rather than a true physical application. 



5. Astrophysical applications in 3D 

To illustrate the possibilities of RAMSES, we present in this 
section two 3D tests of astrophysical significance: the develop- 
ment of Magneto-Rotational Instability (MRI) and the forma- 
tion of a magnetized molecular core. 
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Fig. 7. Time evolution of the magnetic field lines during the current sheet test. The calculation was performed on a uniform 
grid composed of 256 2 cells, using the Roe solver and the MonCen slope limiter. From top left to bottom right, the snapshots 
corresponds successively to times t = 0,0.5, 1, 1.5,2,2.5,3,3.5 and 4. The entire figure can be compared to the figure 12 of 
iGardiner & Stond(l2005al) . 



5.1. The magnetorotational instability 

The development of MHD turbulence resulting from the 
growth of the MRI is likely to be at the origin of an efficient 
radial transport in accretion disks. This instability has been 
extensively studied using finite difference codes like ZEUS 
faawlev et al.lll995h : our goal is here to compare the results 
of RAMSES with those obtained in earlier studies obtained in 
the last decade. 

The MRI is a linear instability that was f irst discovered in 
the 60s dVelikhovl 1 959HChandrasekharl 1 96 ll) before being ap- 
plied to accretion disks theory by Balbus & HawlevNl99lh .lt 
operates in rotating flows threaded by a weak magnetic field 
when the angular velocity decreases outward. Numerical sim- 
ulations applied to accretion disks have shown that the lin- 



ear growth of the instability is followed by MHD turbulence 
that transports angular momentum outward in the disk, thereby 
solving a long standing issue in accretion disk theory (see 
balbus & Hawlevi Il998l for a review). One subclass of these 
simulations has been realised using the so-called shearing box 
approximation. It is a local expansion of the dynamical equa- 
tion in a Cartesian box around a particular radius of the accre- 
tion disk dGoldreich & Lvnden-Belll965h . The interest of this 
local approach is that it is able to capture the dynamics of the 
accretion disk and enables large resolution to be achieved at 
the same time. With such a model, the properties of the MHD 
turbulence can be rather well studied. 

So far, most of the shearing box simulatio ns have been 
done using artificial viscosity codes, like ZEU S jHawlev et alJ 
Il995h . NIRVANA dPapaloizou et alJ |2004 or the Pencil 
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Fig. 8. This figure compares the magnetic field lines obtained at time t = 2 for high-resolution runs (bottom) and low-resolution 
runs (top), as well as for the Roe solver (right) and the local Lax-Friedrich solver (left). Both low and high-resolution runs were 
performed using the AMR scheme: /„„„ = 5 for the former and /„„■„ = 10 for the latter. The refinement strategy is detailed in the 
text. 




Fig. 10. Structure of the flow in the (x,z) plane after 60 orbits. 
The arrows shows the poloidal velocity field overplotted on 
gray scale contours of the y-component of the magnetic field. 
Because of the growth of the MRI, the entire flow has become 
turbulent. 



Code I Brandenburg et ail 1 19951) . Recently. Iciardiner & Stonel 
J2005t ) applied the ATHENA code to the same problem. They 
found that using a Riemann solver make little difference. 

We have implemented the shearing box model in 
RAMSES. To do so, two main modifications have to be made. 
First, the inertial terms appearing in the equations are treated 
as source terms. Second, the boundary conditions need to be 
adapted to the model. They are periodic in y and z, which re- 
spectively correspond to the <p and z coordinates of the stan- 
dard cylindrical coordinates. The boundary conditions in x (the 
equivalent of R in cylindrical coordinates) are the so-called 
shearing boundary conditions widely discussed in the literature 
faawlev et alll995tlGardiner & Stonel2005bl) . 

The initial conditions of our runs are those of the standard 
shearing box model. The initial density is uniform and equal to 
unity everywhere. The velocity is such that 





— qQ,QX 




with q = 1.5 and Qo = 10" 
mal: P = pc^, with cq = 10 



(42) 



3 . The equation of state is isother- 
~ 3 . The initial magnetic field is ini- 
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Fig. 9. This figure shows a zoom on the bottom left magnetic island for the high-resolution run with the Roe solver at time t — 3. 
On the left, only octs boundaries are plotted to clarify the visualization of the AMR grid. On the right, a grey scale image of the 
thermal pressure is shown, with strong transverse shock-waves clearly visible as sharp discontinuities. 
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Fig. 11. Time history of the volume averaged Maxwell stress 
tensor normalized by pressure, obtained in the shearing box 
model. The solid line was computed using RAMSES and the 
dashed line was calculated using ZEUS. Both models use ex- 
actly the same set of parameters. Both shows sustained MHD 
turbulence and a similar amount of angular momentum trans- 
port. 



tially purely vertical and its intensity varies sinusoidally with x 
such that the total net flux threading the computational domain 
vanishes: 



B z = Bq sin 2jtx . 



(43) 



B{) is calculated such that the ratio B between the volume av- 
eraged thermal to magnetic pressure equals 400. The uniform 
grid satisfies -0.5 < x < 0.5, < y < 2iry and -0.5 < z < 0.5. 
The resolution is (N x , N y , N z ) = (32, 100, 32). We ran the model 
with RAMSES, using the Roe Riemann solver and the MonCen 
slope limiter. 

At t = 0, small random velocities are superposed on the 
initial state. The model is evolved during 100 orbital periods 
To, where To = 2n/Q,o- During the first five orbits, the mag- 
netic energy is observed to grow. A measure of the growth 
rate of the instability cr during that period gives cr ~ Q.55Q. 
It is difficult to compare this growth rate with the results of 
linear theory: since the vertical magnetic field varies with x, 
its growth rate does not correspond to that of a single nor- 
mal mode. Nevertheless, we expect the volume averaged evo- 
lution of the magnetic field to be dominated by the growth 
of the field at the position where its strength is initially the 
largest. This is confirmed by a visual inspection of the struc- 
ture of the flow during that phase, which also indicates that 
the associated wavelengths in the radial and vertical direction 
are respectively equal to H a nd H/2. The results of linear the- 
ory dBalbus & Hawlevl 199 lb predicts that cr,/, = 0.55Q in that 
case. Although there is good agreement between cr and cr,/,, we 
want to emphasize that the treatment presented here is only ap- 
proximate. It nevertheless gives confidence in the results of the 
numerical simulation. 

When nonlinear effect becomes important, the magnetic 
energy reaches a peak and start to decline as the whole flow 
breaks down into turbulence, before levelling to a quasi steady 
state it keeps until the end of the simulation. The turbulent na- 
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ture of the disk is illustrated on figure^] It shows the structure 
of the flow in the (x, z) plane after 60 orbits. The arrows rep- 
resent the poloidal velocity field and are overplotted on gray 
scale contours of B Y , which is the dominant component of the 
magnetic field. In order to better quantify the strength of the 
turbulence, we plot on figure ^2 the time history of the vol- 
ume averaged Maxwell stress tensor, normalised by the mean 
thermal pressure Pq: 

rMax < B * B y > 



Po 



(44) 



where < . > denotes a volume average. The solid line in fig- 
ure was obtained with RAMSES. As a comparison, the 
dashed line shows the same quantity obtained with ZEUS for 
the same model (we note that the growth rate measured during 
the linear phase in that case was found to be identical to that ob- 
tained with RAMSES). The two curves are in good agreement 
with each other, even if there is a tendency for RAMSES to 
display some more activity. Indeed, time averaging the curves 
presented on figure^Jbetween 40 and 100 orbits, we obtained 
mean values and rms deviation for the Maxwell stress ten- 
sor equal to (6.6 + 1.5) x 10 3 with RAMSES and equal to 
(3.9 + 1.3) x 10~ 3 with ZEUS. It is worth noting that the (small) 
difference between these two values may not be significant as 
it is an effect of the different dissipative properties of the codes: 
turbulence is driven on large scales by the MRI and damps on 
small scales due to numerical dissipation. The precise saturated 
value of the Maxwell stress results from a balance between 
these two. 

5.2. Magnetized cloud core collapse 

Here we present another 3D test of astrophysical significance: 
the magnetized collapse of a dense prestellar core. In this prob- 
lem the AMR scheme is very useful since the density varies 
over 8 orders of magnitude and the spacial scale, which is about 
the Jeans length, varies over 3-4 orders of magnitude. 

Such calculations in the hydrodynamical case have been 
perfo rmed by several autho rs using either SPH methods 
(e.g. iHe nnebelle et alJ 1200-4 or grid based method (e.g. 
lMatsumoto*&Iianawaf I2003IT In the magnetized case, much 
fewer 3D calculations have been carried out using SPH 
faosking & Whitworfhl 12004. nested grid (Machida et al. 
2005a.b.lBaneriee & Pudritzfeooll) or an AMR implementation 



JZieglerfcOOaT 

In order to do a precise comparison, we adopt the same 
initial conditions as Hosking & Whitworth (2005) and Ziegler 
(2005). The cloud has initially a uniform density of p = 
4.8 x 10~ 18 g cm 4 , a temperature of 10 K and a radius of 
R c = 0.015 pc. The total mass is equal to one solar mass and 
the initial ratio of thermal to gravitational energy is about 0.35. 
The cloud is initially in solid body rotation with angular ve- 
locity a) = 4.25 x 10~ 13 s _1 leading to an initial ratio of ro- 
tational to gravitational energy of about 0.45. To induce frag- 
mentation, an m — 2 perturbation on the density field with an 
amplitude of 10% has been setup initially. The magnetic field 
is initially uniform and parallel to the rotation axis. We use 
the same barotropic equation of state as Hosking & Whitworth 



(2005), namely P = C*pX(l + (p/p crit ) 4/3 ) 1/2 withp crlt = lO" 13 
g cirT 3 . All calculations have initially 64 3 cells and uses 8 addi- 
tional AMR levels. The refinement criteria is based on the Jeans 
length which is numerically described by at least 10 cells. 



5.2.1. Hydrodynamical collapse 

Figure[Q]displays three snapshots of the hydrodynamical case 
performed with the Lax-Friedrich solver. The equatorial den- 
sity field is displayed and the arrows show the velocity field. 
The top panel is very similar to the second panel of Fig. 8 
of Ziegler (2005). Both spiral patterns have approximately the 
same size and the same orientation. Timing is also very similar 
(agreement within about 5% of accuracy). The result compares 
well to the SPH calculation of Hosking & Whitworth (2005) 
shown in their Fig. 2 (bottom-right panel). The second panel of 
Fig. 1 1 21 appears to be similar to the bottom-left panel of Fig. 3 
of Hosking & Whitworth (2005). In both cases about 5 frag- 
ments have formed located approximately at the same position 
in both simulations. The third panel shows the density field 
0.01 Myr after the time of the second panel which is compa- 
rable to the timeshift between left-bottom and right-bottom of 
Fig. 3 of Hosking & Whitworth (2005). The agreement is less 
good than for earlier time. Our results remain symmetric which 
is not the case for the Hosking & Whitworth's results. 

5.2.2. MHD collapse 

In this section we present results for the MHD collapse. The in- 
tensity of the magnetic field is such that the cloud mass-to-flux 
ratio is twice the critical mass-to-flux ratio. This corresponds 
to the first MHD case studied by Ziegler (2005). 

Figure [O] shows results for a calculation performed with 
the Lax-Friedrich solver whereas Fig. [21 shows results ob- 
tained with the Roe solver. Upper panels of Fig.[T3landll4ldis- 
play the equatorial density field whereas bottom panels display 
the density in the x-z plane. Left panels display results at a time 
close to the first panel of Fig. 10 of Ziegler (2005). Right panels 
display results about 0.01 Myr later. 

In his simulation, Ziegler reports the formation of a binary 
having a separation at time t 1 .44% of about 0.06 R c where 
tff is the freefall time. As shown in Fig. [D] there are no bi- 
naries with the Lax-Friedrich solver. There is also a shift in 
time (of about 3%) since the collapse occurs slightly later than 
1.44?//. On the contrary, as shown in Fig. ^] in the simu- 
lation performed with the Roe solver, a binary forms at time 
t = 1.45?// although we find a separation of 0.03 R c instead of 
0.06 R c . Ziegler reports also another case with stronger mag- 
netic field for which the mass-to-flux ratio is 1 .2 times the crit- 
ical mass-to-flux ratio. In this case he finds no binary. We have 
also performed this simulation (not displayed here for concise- 
ness) and we confirm the absence of binary. Although the for- 
mation of the binary appears to be a good numerical test, it 
should be said that it is somehow physically artificial since it re- 
lies on initial conditions for which the density field is perturbed 
but not the magnetic field. As a result, the magnetic Jeans mass 
at the density maximum is lowered, thus making the perturba- 
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tion very unstable. Indeed we have performed a simulation in 
which the m = 2 perturbation has been applied to the magnetic 
field as well. In this case we find no binary. 

Right panels of Fig. H4lshow that the two fragments have 
merged leaving a single central fragment. This is an impor- 
tant difference with the hydrodynamical case in which five 
fragments have been found (although further evolution re- 
veals that two of these fragments merge with the central one). 
Another important departure from the hydrodynamical case is 
the presence of strong outflows (right-bottom panel). There are 
very similar to the outflows shown in Machida et al. (2005b). 
Outflows are also obtained with the Lax-friedrich solver (right- 
botton panel of Fig. I13l> although the flow structure is slightly 
different. With the Roe solver, the velocity field along the z-axis 
vanishes whereas this is not the case with the Lax-Friedrich 



solver. The disk appears to be thicker with the Lax-Friedrich 
solver than with the Roe solver. 

We conclude that RAMSES-MHD is able to reproduce 
quantitatively results obtained by various authors including 
fragmentation and outflows. Significant differences appear be- 
tween results obtained with the Lax-Friedrich and the Roe 
solvers although the former is able to reproduce the main fea- 
tures of the flow. 

6. Conclusion and Perspectives 

In this paper, we have presented an extension of RAMSES to 
MHD. The algorithm is based on the MUSCL-Hancock ap- 
proach already used in the hydrodynamic version of RAMSES 
dTevssierl2002l) . The inducti on equation is evolved in time us- 
ing the standard CT scheme jEvans &Hawle\ll98l . To do so. 
time averaged EMFs are computed on cell edges by solving a 
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2D Ri emann problem, as described in lLondrillo & Del Zannal 
Several tests are presented that illustrate the proper- 
ties and robustness of the code. In particular, we show that the 
AMR scheme implemented in RAMSES can be crucial to de- 
scribe accurately the propagation of some unusual waves pecu- 
liar to MHD like the compound waves. 

We also demonstrate the versatility of RAMSES by study- 
ing two problems of astrophysical significance: the develop- 
ment of MHD turbulence in accretion disk and the collapse of 
dense core in the interstellar medium. In both cases, we report 
results that are consistent with previous studies published in the 
literature. These two applications show that RAMSES is well 
suited to study a wide variety of problems involving MHD in 
astrophysics. 

In future studies, several improvements will now be inves- 
tigated. It will be particularly useful, for example, to develop 
a proper 2D Riemann solver to calculate the time averaged 
EMFs, instead of making linear combination of ID solvers as 



it is done now. Nonlin ear Riemann solvers coul d also be im- 
plemented, like HLLC (IMivoshi & Kusano 2005) for example. 
Obviously, an extension to curvilinear coordinates would also 
be very interesting, particularly for applications involving ac- 
cretion disks or galaxies. Finally, it will be necessary in some 
cases to go beyond the ideal MHD framework and to imple- 
ment new physics like ohmic dissipation or ambipolar diffu- 
sion. 
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Fig. 12. Three timesteps illustrating the hydrodynamical col- 
lapse. 
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